<h2>DESCRIPTION</h2>

<p><em>v.vol.rst</em> interpolates values to a 3-dimensional raster map from
3-dimensional point data (e.g. temperature, rainfall data from climatic
stations, concentrations from drill holes etc.) given in a 3-D vector
point file named <b>input</b>.&nbsp; The size of the output 
3D raster map <b>elevation</b> is given by the current 3D region. Sometimes, the
user
may want to get a 2-D map showing a modelled phenomenon at a
crossection surface. In that case, <b>cross_input</b> and <b>cross_output</b>
options must be specified, with the output 2D raster map <b>cross_output</b>
containing the crossection of the interpolated volume with a surface 
defined by <b>cross_input</b>
2D raster map. As an option, simultaneously with interpolation, 
geometric parameters of the interpolated
phenomenon can be computed (magnitude of gradient, direction of
gradient defined by horizontal and vertical angles), change of gradient,
Gauss-Kronecker curvature, or mean curvature). These geometric
parameteres are saved as
3D raster maps <b>gradient, aspect_horizontal, aspect_vertical, ncurvature, gcurvature, mcurvature</b>,
respectively. Maps <b>aspect_horizontal</b> and <b>aspect_vertical</b> are in degrees.

<p>At first, data points are checked for identical positions and points
that are closer to each other than given <b>dmin</b> are removed.
Parameters <b>wscale</b> and <b>zscale</b> allow the user to re-scale 
the w-values and z-coordinates of the point data (useful e.g. for 
transformation of elevations given in feet to meters, so that the 
proper values of gradient and curvatures can be computed).
Rescaling of z-coordinates (<b>zscale</b>) is also needed when the distances
in vertical direction are much smaller than the horizontal
distances; if that is the case, the value of <b>zscale</b>
should be selected so that the vertical and horizontal distances
have about the same magnitude.

<p>Regularized spline with tension method is used in the interpolation.
The <b>tension</b> parameter controls the distance over which 
each given point influences the resulting volume (with very high tension, 
each point influences only its close neighborhood and the volume goes
rapidly to trend between the points).
Higher values of tension parameter reduce the overshoots that
can appear in volumes with rapid change of gradient. For noisy data, it
is possible to define a global smoothing parameter, <b>smooth</b>.
With the
smoothing parameter set to zero (<b>smooth=0</b>) the resulting volume
passes exactly through the data points. 
When smoothing is used, it is possible to output a vector map <b>deviations</b>
containing deviations of the resulting volume from the given data. 
<p>The user can define a 2D raster map named <b>maskmap</b>, which will
be used as a mask. The interpolation is skipped for 3-dimensional cells
whose 2-dimensional projection has a zero value in the mask. Zero values will
be assigned to these cells in all output 3D raster maps. 
<p>If the number of given points is greater than 700, segmented
processing is used. The region is split into 3-dimensional "box"
segments, each having less than <b>segmax</b> points and interpolation
is performed on each segment of the region. To ensure the smooth
connection of segments, the interpolation function for each segment is
computed using the points in the given segment
and the points in its neighborhood. The minimum number of points taken
for interpolation is controlled by <b>npmin</b> , the value of which
must be larger than <b>segmax</b> and less than 700. This limit of 700 was
selected to ensure the numerical stability and efficiency of the
algorithm. 


<h3>SQL support</h3>

Using the <b>where</b> parameter, the interpolation can be limited to use
only a subset of the input vectors.

<div class="code"><pre>
# preparation as in above example
v.vol.rst elevrand_3d wcol=soilrange elevation=soilrange zscale=100 where="soilrange > 3"
</pre></div>


<h3>Cross validation procedure</h3>

Sometimes it can be difficult to figure out the proper values of
interpolation parameters. In this case, the user can use a
crossvalidation procedure using <b>-c</b> flag (a.k.a. "jack-knife"
method) to find optimal parameters for given data. In this method,
every point in the input point file is temporarily excluded from the
computation and interpolation error for this point location is
computed.  During this procedure no output grid files can be
simultanuously computed.  The procedure for larger datasets may take a
very long time, so it might be worth to use just a sample data
representing the whole dataset.

<p><i>Example 
 (based on <a href="http://www.grassbook.org/data_menu2nd.php">Slovakia3d dataset</a>):</i>
<p><div class="code"><pre>
v.info -c precip3d
g.region n=5530000 s=5275000 w=4186000 e=4631000 res=500 -p
v.vol.rst -c input=precip3d wcolumn=precip zscale=50 segmax=700 cvdev=cvdevmap tension=10
v.db.select cvdevmap
v.univar cvdevmap col=flt1 type=point
</pre></div>

Based on these results, the parameters will have to be optimized. It is
recommended to plot the CV error as curve while modifying
the parameters.
<p>The best approach is to start with <b>tension</b>, <b>smooth</b>
and <b>zscale</b> with rough steps, or to set <b>zscale</b> to a
constant somewhere between 30-60. This helps to find minimal RMSE
values while then finer steps can be used in all parameters. The
reasonable range is <b>tension</b>=10...100,
<b>smooth</b>=0.1...1.0, <b>zscale</b>=10...100.
<p>In <em>v.vol.rst</em> the tension parameter is much more sensitive to
changes than in <em>v.surf.rst</em>, 
therefore the user should always check the
result by visual inspection. Minimizing CV does not always provide the best
result, especially when the density of data are insufficient. Then
the optimal result found by CV is an oversmoothed surface.

<h2>NOTES</h2>
The vector points map must be a 3D vector map (x, y, z as geometry).
The module <a href="v.in.db.html">v.in.db</a> can be used to generate
a 3D vector map from a table containing x,y,z columns.

Also, the input data should be in a projected coordinate system, such as
Universal Transverse Mercator. The module does not appear to have support for 
geographic (Lat/Long) coordinates as of May 2009.

<p><em>v.vol.rst</em> uses regularized spline with tension for
interpolation from point data (as described in Mitasova and Mitas,
1993). The implementation has an improved segmentation procedure based
on Oct-trees which enhances the efficiency for large data sets. 

<p>Geometric parameters - magnitude of gradient (<b>gradient</b>),
horizontal (<b>aspect_horizontal</b>) and vertical (<b>aspect_vertical)</b>aspects,
change of gradient (<b>ncurvature</b>), Gauss-Kronecker (<b>gcurvature</b>) and
mean curvatures (<b>mcurvature</b>) are computed directly from the
interpolation function so that the important relationships between
these parameters are preserved. More information on these parameters
can be found in Mitasova et al., 1995 or Thorpe, 1979.

<p>The program gives warning when significant overshoots appear and
higher tension should be used. However, with tension too high the
resulting volume will have local maximum in each given point
and everywhere else the volume goes rapidly to trend. With a smoothing
parameter greater than zero, the volume will not pass through the data 
points and the higher the parameter the closer the volume will be to the 
trend. For theory on smoothing with splines see Talmi and Gilat, 1977 or Wahba, 1990. 

<p>If a visible connection of segments appears, the program should be
rerun with higher <b>npmin</b> to get more points from the
neighborhood of given segment. 

<p>If the number of points in a vector map is less than 400, <b>segmax</b>
should be set to 400 so that segmentation is not performed when it is
not necessary. 

<p>The program gives a warning when the user wants to interpolate outside the
"box" given by minimum and maximum coordinates in the input vector map. 
To remedy this, zoom into the area encompassing the input vector data points.

<p>For large data sets (thousands of data points), it is suggested to
zoom into a smaller representative area and test whether the parameters
chosen (e.g. defaults) are appropriate. 

<p>The user must run <em>g.region</em> before the program to set the
3D region for interpolation. 


<h2>EXAMPLES</h2>

<!-- TODO: find better data. This example is nonsensical :-) -->
Spearfish example (we first simulate 3D soil range data):

<div class="code"><pre>
g.region -dp
# define volume
g.region res=100 tbres=100 res3=100 b=0 t=1500 -ap3

### First part: generate synthetic 3D data (true 3D soil data preferred)
# generate random positions from elevation map (2D)
r.random elevation.10m vector_output=elevrand n=200

# generate synthetic values
v.db.addcolumn elevrand col="x double precision, y double precision"
v.to.db elevrand option=coor col=x,y
v.db.select elevrand

# create new 3D map
v.in.db elevrand out=elevrand_3d x=x y=y z=value key=cat
v.info -c elevrand_3d
v.info -t elevrand_3d

# remove the now superfluous 'x', 'y' and 'value' (z) columns
v.db.dropcolumn elevrand_3d col=x
v.db.dropcolumn elevrand_3d col=y
v.db.dropcolumn elevrand_3d col=value

# add attribute to have data available for 3D interpolation
# (Soil range types taken from the USDA Soil Survey)
d.mon wx0
d.rast soils.range
d.vect elevrand_3d
v.db.addcolumn elevrand_3d col="soilrange integer"
v.what.rast elevrand_3d col=soilrange rast=soils.range

# fix 0 (no data in raster map) to NULL:
v.db.update elevrand_3d col=soilrange value=NULL where="soilrange=0"
v.db.select elevrand_3d

# optionally: check 3D points in Paraview
v.out.vtk input=elevrand_3d output=elevrand_3d.vtk type=point dp=2
paraview --data=elevrand_3d.vtk

### Second part: 3D interpolation from 3D point data
# interpolate volume to "soilrange" voxel map
v.vol.rst input=elevrand_3d wcol=soilrange elevation=soilrange zscale=100

# visualize I: in GRASS GIS wxGUI
g.gui
# load: 2D raster map: elevation.10m
#       3D raster map: soilrange

# visualize II: export to Paraview
r.mapcalc "bottom = 0.0"
r3.out.vtk -s input=soilrange top=elevation.10m bottom=bottom dp=2 output=volume.vtk
paraview --data=volume.vtk
</pre></div>

<h2>KNOWN ISSUES</h2>
<b>deviations</b> file is written as 2D and deviations are not written as attributes.

<h2>REFERENCES</h2>
<p>Hofierka J., Parajka J., Mitasova H., Mitas L., 2002, Multivariate
Interpolation of Precipitation Using Regularized Spline with Tension.
Transactions in
GIS&nbsp; 6, pp. 135-150. 
<p><a href="http://www4.ncsu.edu/~hmitaso/gmslab/">Mitas, L.,
Mitasova, H.</a>, 1999, Spatial Interpolation. In: P.Longley, M.F.
Goodchild, D.J. Maguire, D.W.Rhind (Eds.), Geographical Information
Systems: Principles, Techniques, Management and Applications, Wiley,
pp.481-492 
<p>Mitas L., Brown W. M., Mitasova H., 1997,
<a href="http://www4.ncsu.edu/~hmitaso/gmslab/lcgfin/cg-mitas.html">Role
of dynamic cartography in simulations of landscape processes based on
multi-variate fields.</a> Computers and Geosciences, Vol. 23, No. 4,
pp. 437-446 (includes CDROM and WWW: www.elsevier.nl/locate/cgvis) 
<p>Mitasova H., Mitas L.,&nbsp; Brown W.M.,&nbsp; D.P. Gerdes, I.
Kosinovsky, Baker, T.1995, Modeling spatially and temporally
distributed phenomena:
New methods and tools for GRASS GIS. International Journal of GIS, 9
(4),
special issue on Integrating GIS and Environmental modeling, 433-446. 
<p> Mitasova, H., Mitas, L., Brown, B., Kosinovsky, I., Baker, T.,
Gerdes, D. (1994):
<a href="http://www4.ncsu.edu/~hmitaso/gmslab/viz/ches.html">Multidimensional
interpolation and visualization in GRASS GIS</a> 
<p><a href="http://www4.ncsu.edu/~hmitaso/gmslab/papers/lmg.rev1.ps">Mitasova
H. and Mitas L. 1993</a>: Interpolation by Regularized Spline with
Tension: I. Theory and Implementation, <i>Mathematical Geology</i> 25,
641-655. 
<p><a href="http://www4.ncsu.edu/~hmitaso/gmslab/papers/hmg.rev1.ps">Mitasova
H. and Hofierka J. 1993</a>: Interpolation by Regularized Spline with
Tension: II. Application to Terrain Modeling and Surface Geometry
Analysis, <i>Mathematical Geology</i> 25, 657-667. 
<p>Mitasova, H., 1992 : New capabilities for interpolation and
topographic analysis in GRASS, GRASSclippings 6, No.2 (summer), p.13. 
<p>Wahba, G., 1990 : Spline Models for Observational Data, CNMS-NSF
Regional Conference series in applied mathematics, 59, SIAM,
Philadelphia, Pennsylvania. 
<p>Mitas, L., Mitasova H., 1988 : General variational approach to the
interpolation problem, Computers and Mathematics with Applications 16,
p. 983 
<p>Talmi, A. and Gilat, G., 1977 : Method for Smooth Approximation of
Data, Journal of Computational Physics, 23, p.93-123. 
<p>Thorpe, J. A. (1979): Elementary Topics in Differential Geometry.
Springer-Verlag, New York, pp. 6-94.

<h2>SEE ALSO</h2>

<em>
<a href="g.region.html">g.region</a>,
<a href="v.in.ascii.html">v.in.ascii</a>,
<a href="r3.mask.html">r3.mask</a>,
<a href="v.in.db.html">v.in.db</a>,
<a href="v.surf.rst.html">v.surf.rst</a>,
<a href="v.univar.html">v.univar</a>
</em>

<h2>AUTHOR</h2>

Original version of program (in FORTRAN) and GRASS enhancements: <br>
Lubos Mitas, NCSA, University of Illinois at Urbana-Champaign,
Illinois, USA, since 2000 at Department of Physics, 
North Carolina State University, Raleigh, USA
<a href="mailto:lubos_mitas@ncsu.edu">lubos_mitas@ncsu.edu</a>
<br>
Helena Mitasova, Department of Marine, Earth and Atmospheric Sciences,
North Carolina State University, Raleigh, USA,
<a href="mailto:hmitaso@unity.ncsu.edu">hmitaso@unity.ncsu.edu</a>
<p>
Modified program (translated to C, adapted for GRASS, new
segmentation procedure): <br>
Irina Kosinovsky, US Army CERL, Champaign, Illinois, USA <br>
Dave Gerdes, US Army CERL, Champaign, Illinois, USA 
<p>
Modifications for g3d library, geometric parameters,
cross-validation, deviations: <br>
Jaro Hofierka, Department of Geography and Regional Development,
University of Presov, Presov, Slovakia,
<a href="mailto:hofierka@fhpv.unipo.sk">hofierka@fhpv.unipo.sk</a>,
<a href="http://www.geomodel.sk">http://www.geomodel.sk</a>

<!--
<p>
<i>Last changed: $Date$</i>
--> 
